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A Computer Subroutine for Stress 
Analysis of Rotating, Heated Disks 
by 

John E. Brock and Robert E. Brown 
Introduction 

Although, as is indicated by the title hereof, the principal 
purpose of this monograph is to present tested and proved digital com- 
puter software for the analysis of stress in a spinning axisymmetric 
disk having a radially variable thermal strain field, the opportunity 
is also taken of developing the theory and presenting some analytic 
solutions . 

The method developed herein for computer analysis of disks 
having a general law of thickness variation was suggested by the al- 
gorithm contained in reference 2 and it appears to have advantages 
over such procedures as that of M. Donath, reference 3, which has been 
widely circulated in a book by S. Timoshenko, reference 5. 

In what follows we immediately obtain a second order linear 
differential equation with dependent variable u, the radial deformation, 
and r, the radius. Analytic treatment is given for two particular laws 
of thickness variation. For the general case of thickness variation, 
the equation is recast as a second order linear differential equation 
in which the dependent variable is radial stress, c? r . However, for 
numerical treatment an alternate form is more useful and direct and 
this forms the basis of the digital computer software ’which is given 



and illustrated in the appendices hereto. 



Fundamental Analysis 

We presume that the disk is thin enough and that the thick- 
ness varies slowly enough with respect to radius that we are just- 
ified in neglecting all stress components excepting only the radial 
stress o r and the circumferential stress cr 0 . Material properties 
E, Young's modulus, and v, Poisson's ratio, are presumed to be in- 
deed constant. The thermal strain field, aT, and the density y 
may be specified functions of radius. 

Two types of problem are considered: 

1. Annular disk, 0 < a £ r = b, with a (a) and o (b) 

r r 

being specified. 

2. Solid disk, 0 i r = b, with c^Cb) being specified. 

We also use the symbols t = t(r) for disk thickness and w for 
angular velocity. Other simplifying notation will be introduced 
later on. 

Consideration of radial equilibrium leads without difficulty 
to the equation 



1 d(rta ) 
dr 



a Q - yoj 2 r 2 



( 1 ) 



The thermoelastic constitutive equations are 

Ee. = a Q - va + EaT; Ee = a - vcn + EaT (2a, b. 

y y v v v t? 

where the strain components are 

e Q = u/r; e r = du/dr (3a,b 

Strain cornpatability leads to the equation 

r ^(Eu/r) = -(l+v)(a 0 -a r ) (4) 

These equations may easily be combined into the differential equation 
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u" + u'/r - u/r 2 + (t'/t)[u' + vu/r - (l+v)ctT] 

= (l+v)(aT) ' - (1-v 2 )yco 2 r/E (5) 

where primes denote differentiation with respect to r. 

Two particular laws of thickness variation permit simple 
analytic treatment. 

Power Law of Thickness Variation 
If 

t = t 0 (r/r 0 ) n (6) 

so that 

(t'/t) = n/r (7) 

then equation 5 becomes 

u" + (l+vn)u'/r + (vn-l)u/r 2 = 8' + nB/r - kr (8) 

where 

8 = (l+v)cxT, k = (l-v 2 )yw 2 /E (9,10) 

Equation 8 may be rewritten as 

d r 1-m d , (m+n)/2 (2+n-m)/2 

gp[r 5r^ r u) ] = r (8' + n8/r -kr) (11) 

where 

m = ±/(n 2 -4 n+4) = ±/[(n-2) 2 +4(l- )n] (12) 

and either the positive or the negative sign may be used Equation 11 
may be proved simply by performing the indicated operations and com- 
paring with equation 8. 

The quantity on the right in equation 11 is well defined so 
that the solution of the differential equation may be obtained s irmly 
by integration, multiplication by r m- ‘*', and another integration. Two 
constants of integration are introduced. For the solid disk (case 2) 
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u(0) = 0 gives one condition and the second comes from satisfying 
the given value of a r (b). For the annular disk (case 1) satisfying 
the given conditions a^(a) and a^(b) permits evaluating the con- 
stants . 

An example with n = - 0.42 is given in Appendix C. Note 
that the case n = 0 corresponds to a disk of uniform thickness. 



Exponential Law of Thickness Variation 



If 

t = t Q exp(-mr 2 ) 

where m is a constant of appropriate dimensionality, then 
(t'/t) = -2rm 
and equation 5 becomes 

u" + (1/r - 2rm)u' - (1/r 2 + 2vm) = -2rm8 + 8' - kr 
If additionally we assume that 8' = 0 (which makes the thermal 
strain field constant — a triviality) the solution is simply 



(13) 

(14) 

(15) 



u = r(B+k/2m)/(l+v) 



( 16 ) 



In this case we cam easily find 

a r = a Q = ya> 2 /2m (17) 

which is independent of r. Thus, if an allowable normal stress a A 
is specified and if blade or bucket loading at r = b is w (oounds, 
say) per unit circumference, then a disk having thickness 

t = (w/o A ) exp[(b 2 - r 2 )(yw z /2a A )] ( 18 ) 

will be such that a r = a Q = a A « If the failure criterion is the 
maximum shearing stress criterion (Tresca's condition), it is clear 
that this disk is optimal in the sense of having least volume and 



thus having least weight. If the failure criterion were that of 



von Jiises, a slightly lighter disk would suffice. 

General Law of Thickness Variation 
Equations 3a and 2a give 

Eu = rEaT + r(ag-vo r ) (19) 

and by differentiation 

Eu' = EaT + rE(aT) ' + (a„-va ) + r (a ’-vo ') (20) 

t) r o r 

From 3b and 2b we also have 

Eu' = EaT + a - va Q (21) 

r t) 

Subtracting 21 from 20 and rearranging gives 

E(aT) ' + (l+v)(a Q -a r )/r + a Q ' - va r ' = 0 (22) 

Equation 1 may be rewritten as 

a fi -a = ra ' + rvo + yoj 2 r 2 (23) 

t) V V V 

where we have written 

v = t'/t (2*1) 

for convenience. Differentiating 23 we get 

o q = 2o r ' + ra r " + va r + rvo r ' + rv'a r + 2yw 2 r (25) 

and substituting 23 and 25 into 22 gives 



ro r " + (3wv)a r ' + f(2tv)vm']c r + C3+v) Y » ! r + EM)- - 0 (26) 

This is a single differential equation with dependent var- 
iable a and can be dealt with by standard numerical methods. The 
r J 

conditions for evaluating the constants of integration have been 
mentioned earlier. 

However, the preceding procedure is not particularly satis- 
factory. For one thing, the solid disk case for which r can become 
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zero encounters numerical difficulties unless special treatment is 
employed to evade them. More importantly, however, if T is given 
by a graph or a numerical table, determination of (aT) ' may involve 
numerical difficulties which, ultimately, are due to our having per- 
formed the differentiation to arrive at equation 20. Accordingly, 
an alternate procedure which adheres more closely to the fundamental 
mechanics of the problem is described in what follows. 

We consider the annular disk first and we represent the un- 
known stress difference a. -a in the form 

0 r 

a 0^r = (A+Bn)r (27) 

where A and B are unknown constants and n is an unknown function 
normalized so that 



n(a) = 0, n(b) = 1 (28a, b) 

Initially we make an assumption for n taking -a linear var- 
iation in the absence of better information. By the use of some of 
the preceding equations we will be able to construct an improved 
form for n and will iterate until there is satisfactory convergence. 



Let 



z = Eu/r, w = EaT, B = yco 2 (29,30,31) 

noting that w and 3 have different meanings here than when they 

were used earlier. We can recast equation 1 as 

d(ta r )/dr = (a Q -o r )t/r - yco 2 rt = t(A+Bn-yaj 2 r) (32) 

Before performing the indicated integration we introduce two 

convenient notational devices, viz. 

r r 



p. . . 




c • • • ] r=b 



(33,34) 



Thus, from equation 32 we have 
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( 35 ) 



ta = (to ) + Apt + Bppt - p8rt 

r* r 3, 

From equation 4 and equation 27 we have 

dz/dr = -(1+v) (A+Bn ) (36) 

so that 

z = (z) - (l+v)[A(r-a) + Bpn] (37) 

a. 

Equation 2a is 

z = w + (a 0 — o r ) + (l-v)a r = w + (A+Bn)r + (l-v)a^ (38) 

so that 



(z) = (w) 


+ Aa + (l-v)(a ) 


(39) 


a a 


r a 


(Z) b = <w) b 


+ (A+B)b + (l-v)(a r ) b 


(40) 



Evaluating 38 at r = b and using 39 and 40 gives one equation 
involving the unknowns A and B. Evaluating 35 at r = b gives a 
second. These equations can be arranged as 



*(pt) *(pr|t) 


A' 


= 


'*(pSrt) + (ta r ) b - (ta r ) a 


(2+v)(b-a) b+(l+v)*(pn) 


_B_ 




.(w) a -(w) b+ (l-v)[(o r ) a -(<, r ) b ]. 



so that one can easily solve for A and B. Then is obtained from 

35, z is obtained from 37 and 39, and (a Q -a r ) and a g are obtained 

from 38. Then a new function n is calculated from 

n = [(a Q -o r )/r-A]/B (42) 

Using the new n the process is iterated, convergence being 

monitored by examination of the sequence of values of A and B that 

are calculated. When convergence is satisfactory, the desired 

functions a and a Q are at hand, 
r y 

The situation is simpler for the solid disk, case 2. (o r ) a 
is not given but conditions of continuity require that (a g -a r )/r 
vanish at r = 0. Thus A is zero and B can be obtained from 
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( 43 ) 



(w) 0 - (w) b + c[l-(t)o/(t) b ](ta r ) b + c*(p3rt) 

B = b +' (1+V)*rpn7 + ' c^Cpht J 

where 

c = (l-v)/(t)o (44) 

Otherwise the procedure is as for case 1. 

Computer Implementation 

The theory embodied in equations 27 through 44 and the 
associated procedure has been programmed for digital computer 
using the FORTRAN language. An initial programming based directly 
on the preceding equations and written in January 1978 by the 
junior author hereof has been supplanted by a newer programming 
which is somewhat more compact due to the employment therein of 
ancillary subroutines developed for use in another problem (the 
lateral buckling of elastic beams) on which we are working. This 
program, actually a subroutine called RODISK is listed in Appen- 
dix A hereof. This listing itself contains comments which ade- 
quately explain the construction of a MAIN program which supplies 
necessary input information and which invokes RODISK. Appendix 
B lists the ancillary subroutines. In each case comments indicate 
the purpose and employment of the subroutine. These may prove 
useful in constructing the MAIN program. For this reason, the 
ancillary subroutines DUFV and PRIV are given even through they 
are not called by RODISK. 

In the computer implementation, the various functions of 
r which appear in the theory are represented by vectors the ele- 
ments of which are function values at equally spaced values of r. 
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The ancillary subroutines manipulate these vectors. All of these 
are obvious except possibly INTV which performs an integration by 
use of Milne's formulas, cf. reference 4. 

In the subroutine RODISK there is a slight departure from 
the theory as given herein. As a first step, all quantities and 
functions were "dedimensionalized" but otherwise the method is 
just as described above. Somewhat finer details of what ’was 
actually done may be found in reference 1. 

Appendix C contains some examples and remarks concerning 

them. 
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Appendix A 



Listing of subroutine RODISK 

(The listing on this page, page 10, is of the comments which 

provide instructions for the use of RODISK. Commands appear 

on the following two pages . ) 

SUBROUTINE RODISK. JOHN E. BRGCK., 1 MAY 1978 
THIS IS A SUBROUTINE FOR DETERMINING RADIAL AND CIRCUM- 
FERENTIAL STRESSES IN AN AXISYMMETRIC THIN ELASTIC DISK, 
ROTATING AT ANGULAR VELOCITY CM5GA (RADIANS PER SECOND) 
ABOUT ITS AXIS OF SYMMETRY AND HAVING AN AXISYMMETRIC 
DISTRIBUTION CF THERMAL STRAIN. TWO TYPES OF PROBLEM MAY 
B'E TREATED: 

TYPE 1. ANNULAR DISK OF INSIDE RADIUS ARAD AND 
OUTSIDF RAD T US 8RAD. THE RADIAL STRESS IS S RA AT THE 
INNER RADIUS AND SRB AT THE OUTER RADIUS. THE INSIDE 
RADIUS MUST BE GREATER THAN ZERO. 

TYPE 2. SOLID DISK OF OUTSIDE RADIUS BRAD AND WITH 
RACIAL STRESS SRB AT THE OUTSIDE RADIUS. 

THE USER MUST PROVIDE A MAIN PROGRAM WHICH CALLS SUBROUTIN 
RODISK AFTER IT HAS SUPPLIED THE FOLLOWING INFORMATION. 

(L> N, INTEGER. C N— 1 ) IS THE NUMBER OF EQUAL SUBDIVISION 
INTO WHICH THE ANNULAR RACIUS (5RAC MINUS ARAD) IS 
DIVIDED FOR COMPUTATIONAL PURPOSES. T HE PRESENT 
DIMENSIONING CAN ACCOMMODATE N NOT LARGER THAN 101. 

(2) BRAD 

(3) ARAD (NOT NECESSARY FOR PROBLEMS CF TYPE 2.) 

(4) SRB 

(5) SR A (NOT NECESSARY FOR PROBLEMS OF TYPE 2.) 

(6) TEEBEE, DISK ’’HICKNESS AT OUTSIDE RADIUS 

(7) POIS, POISSON'S RATIO. 

(8) KP ( 1 ) = 1,2. INTEGER TO DENOTE PROBLEM CF TYPE 1,2. 

(9) KP ( 2 ) , INTEGER TO PROVIDE FOP SKIPPING WHILE PRINTING 
OUTPUT. FOR EXAMPLE, IF N = 101 AND KP(2)=5, ONLY 
EVEPY FIFTH SET OF VALUES WILL BE PRINTED: 1ST, 6TH , 
..., 96 T H, 1 01 ST . 

(10) KP ( 3) , INTEGER SPECIFYING the NUMBER OF ITERATIONS 
TO BE PERFORMED. USUALLY KP(3)=10 IS SUFFICIENT FOR 
ENGINEERING ACCURACY. 

(11) KP ( 4 ) . IF K P ( 4) =0 ONLY FINAL ANSWERS WILL eE PRINTED 

IF KP ( 4 ) = 1 A SEQUENCE OF ITERANT VALUES WILL BE 
PRINTED, INDICATING DEGREE OF CONVERGENCE. 

(12) VECTORS X < I , J ) , 1 = 1,2, 3; J=1,2,...,N. 

VECTOR X(1,J) CONTAINS VALUES 0>= THE RATIO (LOCAL 
THICKNESS n F D I SK ) / ( T EE BE E > COMPUTED AT EQUALLY SPACED 
RADII STARTING AT THE INSIDE AND ENDING AT THE OUTSIDE. 

VECTOR X ( 2 , J ) CONTAINS VALUES OF ( GAMM A ) ( OMEGA-SQUA R E 
TIMES ( BRAD- SCAR ED ) DIVIDED BY (SRB). FOR MOST PROBLEMS 
GAMMA DOES NOT VARY WITH RADIUS AND THIS QUANTITY IS A 
CONSTANT. 

VECTOR X ( 3 , J ) CONTAINS VALUES OF ( F )( ALPHA ) (TE E )/ (SRB 
WHERE (E) IS YOUNG'S MODULUS, (ALPHA) IS THE COEFFICIENT 0 
LINEAR THERMAL EXPANSION, AND (TEE) IS TEMPERATURE CHANGE. 

THE MAIN PROGRAM MUST CONTAIN TH C STATEMENTS: 

IMPLICIT REAL*8 (A-H,0-Z) 

INTEGER KP ( 4 ) 

COMMON X ,N 

COMMON /ONE/ A RAD, BR AD , S RA , SR B , T EEBEE , PO IS,KP 

FOLLOWING SUBROUTINE RODISK THERE ARE SEVERAL 
ANCILLARY SUBROUTINES WHICH PERFORM VARIOUS OPERATIONS 
ON THE VECTOFS X(I,J). THE FUNCTION OF EACH IS OBVIOUS 
FROM THE LISTING. THEY MAY BE USED IN THE USER'S MAIN 
PROGRAM. TWO CF THESE ANCILLARY SUBROUTINES WHICH ARE 
AVAILABLE IN THIS PACKAGE BLT WHICH ARE NOT CALLEC BY 
SUBROUTINE PODISK ARE DUPV WHICH DUPLICATES A VECTOR AND 
PR I V WHICH PRINTS A VECTOR. 
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SUBROUTINE ROD I SK 
IMPLICIT R2AL*3 (A-H, O-Z ) 

REALMS X( 20 ,1 01 ) 

INTEGER KP ( 4 ) 

COMMON /ONE /ARAD, BRAD ,SRA ,SRB»T 
ON E = I .D+0 
PO 1 3 = 2 . D“ 1 
Z p RO = C .0+0 

I F (KP ( 1 ) . EQ .2 ) ARAD = ZERO 
RHC = ARAD/ BRAD 



E8EE,PDIS,KP 



ENM=N-1 , , 

WRITE (6 ,2 ) KP ( 1 ) 

2 FORMAT ( // , I OX , ' ROD I SK PROBLEM OF TYPE ',11,' 
DG 5 1=1, N 
E I M= I - 1 
Y = EI M/5NM 

X ( 4 , I ) =RHO+ ( ONE -PHO ) *Y 
X (5,1 )=Y 
5 X (6,1 )=Y 
I TER = I 

ETA = X(1,1)/X(1,N) 

IF(KP (1 ). EQ.2) GO TO 100 

THE PROBLEM IS OF TYPE ONE : ANNULAR DISK 
SPAT = SRA/SR B 

C 1 = ( 2 . D+0+ POI S ) * ( ONE-RHO ) 

CALL INTV (1,7) 

C2=X ( 7,N) 

C5=X(3,1)-X(3,N)-(ONE-P0IS)*(ON=-SRAT ) 

CALL MUL V ( 1,2,8) 

CALL MU LV ( 8 ,4 , 9 ) 

CALL I NTV ( 9 , 1 0 ) 

C6=X(1C,N)+(0NF-ETA*SRAT)/(0NE-RH0 ) 

20 CALL INTV(6 ,11 ) 

C3=0NE+(0NZ-RH0)*(0NE+P0IS)*X(1 1,N) 

CALL MULV ( 1,6,12) 

CALL I NTV ( 1 2 , 1 3 ) 

C 4 = X ( 13, N) 

D = C1*C4-C2*C3 
A = (C5*C4-C6*C3 ) /D 
B= (C 1 *C6-C 2*C 5 ) /O 

IF(KP(4).EQ.l) WRI T E ( 6 , 7 ) ITER, A, 3 
7 F0RMAT(5X,I 10 ,1P2E20.5) 

CALL MULS (7 , 14 , A ) 

CALL MU LS (13, 15 ,3) 

CALL ADDV< 14, 15,15) 

CALL S'JBV (15,10,16) 

S=CNE-RHO 

CALL MULS ( 1 6 , 1 6 , S ) 

S = ET A*SRAT 

CALL ADDS (16, 16, S) 

CALL DI VV( 16, 1, 20) 

ZJ=A*RHO+X( 3, 1 ) +(ONE-PCIS )*SRAT 
CALL MULS ( 11 ,11 ,3) 

CALL MULS ( 5 , 16 , A) 

CALL ADDV ( 1 1, 16 ,16) 

S =- l ONE-RHO ) * ( ONE+POI $) 

CALL MULS( 16,16, S) 

CALL ADDS ( 1 6 , 19 , ZO ) 

S=POI S-ONE 

CALL MULS (20,18,5) 

CALL ADDV( 18,19,18) 

CALL SUBV( 18,3, 19) 

I TER = ITER + 1 

I F (I TER.GT. KP (3 ) ) GO TO 200 
CALL DIVV( 19,4, 17) 

CALL SUBS ( 17, 17 ,A ) 

CALL CIVS(17,6,B) 

GO TO 20 



•,//) 
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ioo ci=pcis-cn: 

C 2=ET A 

C5 =X ( 3 , 1 ) -X ( 3 » N ) +PO I S -CNF: 

CALL NU LV ( 1 » 2 , 7 ) 

CALL VULV (7,4,8) 

CALL INTV (6,9) 

C 6 =3NE + X ( 9 » N ) 

105 CALL INTV (6,10) 

C3=ONE+(ONE+POIS)*X( 10,N) 

CALL NULV( 1 1 6 * 1 1 ) 

CALL INTV (11,12) 

C4=X( 12, N) 

D = C 1 ~ C 4-C 2 * C 3 
S R AT = (C5*C4-C6*C3 )/D 
0 = (C1*C6-C2*C5 ) /D 

I F ( KP ( 4 ) • EQ . 1 ) kRI TE ( 6 ,7 ) ITER, SRA T , S 
CALL VULS ( 12, 13,9) 

CALL SUBV ( 1 3 ,9 , 14) 

S =ET A*SRAT 

CALL A00S ( 14 , 20 , S ) 

CALL DIVV (20,1,20) 

S=-3*(GNE+P0IS) 

CALL PULS( 10,15, S) 

S = P Q I S-ONE 

CALL NULS(20,18,S> 

S = X ( 3 » 1)+(0NE-P0IS) *SR,AT 
CALL A DOS ( 1 8 , 19 , S ) 

CALL SU9V( 19,3,19) 

S=-B* (ONE+POIS ) 

CALL f / ULS(10,17»SI 
CALL ADDV< 19, 17 ,19) 

ITER=ITER + 1 

I F (I T5R.GT . KP ( 3 ) ) GO TC 200 
CALL MUL S (4,16,8) 

X (16 , 1 ) =ONE 

CALL Cl VV (19, 16, 6) 

GO TO 105 

200 CALL NUL$(4 ,7 ,6RAD) 

CALL M ULS ( 20 , 8 , SR6 ) 

CALL MULS (19,9, SR8) 

CALL ADDV (8,9,9) 

S=SRB/8RAD**2 
C ALL RULS ( 2 » 10 » S ) 

CALL PULS<3 ,U,SRB» 

CALL MULS( 1,15, TEEBcE ) 

W P I T E (6,204 > 

204 FORMAT (////) 

WRITE(6,205 > 

235 F ORM AT ( 24 X, 'RADIUS' ,1 IX, 'THICKNESS ' ,6X, • GAMM A . OMEGA . SQ 
1' EE. ALPHA. TEE' »8X,' SIGMA. p ADIAL' ,7X, *S IGMA.C IRCUMF' ) 
KSKI P=KP (2 ) 

DC 210 I=1»N,KSKIP 
J=I/KSKI P 

WRITE(6,211)J,X(7,I),X(15,I),X(1),I) ,X(11,I) ,X( 6,1) * X < 9 

210 CONTINUE 

211 FORMATt 1 10, 1P6E20. 5 ) 

RETURN 

END 
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Appendix B 



Listing of ancillary subroutines 



SUBROUTINE A0DV(N1 , N2 , N3 ) 

ANCILLARY SUBROUTINE 

ADC TWO VECTORS TERM 3Y TERM 

X ( N3 * I ) = X ( N1 ,1 ) +X ( N2 » I ) 

REAL *8 X ( 2 D , 1 0 1 )»S 
COMMON X » N 
DO 1 1*1, N 

I X ( N3 , I ) =X ( N 1 , 1 ) +X ( N2 , I ) 

RETURN 

END 

SUBROUTINE SUBV (N1,N2 ,N3) 

ANCILLARY SUBROUTINE 

SUBTRACT TWO VECTORS T ERM BY term 

X ( N 3 , I ) = X( N 1 , 1 ) - X ( N 2 , 1 ) 

REALMS X<2J,101),S 
COMMON X, N 
DO 1 1 = 1, N 

1 X ( N3 , I ) = X ( N 1 , I ) -X C N 2 , I ) 

RETURN 

END 



SUBROUTINE MU L V ( N 1 , N 2 , N 3 ) 

ANCILLARY SUBROUTINE 

MULTIPLY TWO VECTORS T ERM 3Y TERM 

X(N3,I)=X(N1,I ) *X ( N2 , I ) 

REAL* 8 X ( 2 0 » 1 0 1 ) , S 
COMMON X» N 
DC 1 1=1, N 

X ( N3 , I )=X(N1, I )*X(N2 , I ) 

RETURN 

END 



1 



SUBROUTINE DIVV(N1,N2,N3) 
ANCILLARY SUBROUTINE 
DIVIDE TWO VECTOR :> iFkM dY 
X (N3 ,1 ) =X(N1 , I )/X( N2, I ) 
IEAL* £ X( 20 , 101 ) , 5 
COMMON X , N 
DC 1 1=1, N 

< l N 3 . I ) =X ( N 1 , 1 ) /X ( N2 , I ) 



RETURN 

ENC 



TERM 



1 



SUBROUTINE ADDS ( N1 . N2 » S ) 
ANCILLARY SUBROUTINE 
ADC A SCALAR TO EACH 
X ( N2 , I ) = X ( Ni ,1 H-5 
REALMS X( 20 , 1 01 ) , S 
COMMON X,N 
DO 1 1=1, N 
X ( N2 , I ) = X(N1. I)+S 



RETURN 

END 



TERM OF A VECTOR 



- 13 - 



oo ooo ooo ooo ooo 



SUBROUTINE SUBS ( N 1 •» N 2 , S ) 

ANCILLARY SUEROUTINt 

SUBTRACT A SCALAR FROM EACH TERM GF A VECTOR 
X (N2,I >=X(N1 ,1 )-S 
REALMS X(20,L0L),S 
COMMON X , N 
00 1 I = 1 , N 
1 X (N2, I )=X(N1,I ) -S 
RETURN 
END 

SUBROUTINE MU L S ( N1 , N2 , S ) 

ANCILLARY SUBROUTINE 

MULTIPLY EACF TERM OF A VECTOR 8Y A SCALAR 
X ( N2 » I ) - X( N1 ,1 )*S 
R E AL* 8 X( 20 , 1 01 ) » S 
COMMON X * N 
00 1 1=1, N 
1 X ( N2 , I ) = X { N 1, 1 > *S 
RETURN 
END 

SLBROLTINE DI VS (N1 , N2 ,5 ) 

ANCILLARY SUBROUTINE 

DIVIDE EACH TERM OF A VECTOR BY A SCALAR 
X ( N 2 , 1 ) = X( N1 , 1 J / S 
REAL*8 X( 20 , 10 1 ) , S 
COMMON X , N 
DG 1 1=1, N 
1 X (N2, I ) = X ( N 1 , I ) /S 
RETURN 
END 

SUBROUTINE 0UPV(N1,N2) 

ANCILLARY SUBROUTINE 
DUPLICATE A VFCTOR 
X (N2,I ) = X ( N 1 , 1 ) 

REAL* 8 X( 20 , 1 01 ),S 
COMMON X , N 
DO 1 1=1, N 
1 X (N2 , I ) =X (N 1 , I ) 

RETURN 

END 

SUBROUTINE INTV(N1,N2) 

ANCILLARY SUBROUTINF 

INTEGRATE A VECTOR USING MILNE'S METHOD 
REAL* 8 X(20 ,101 ) , EN , R , ADD , N I NO, NTNO , F I VO ,T HTO 
COMMON X , N 
E N=N- 1 
EN = 1 . O + O/EN 
R=EN/2.4D+1 
N I N0=R*9.D+0 
NTN0=R*1.9D+1 
F I VQ = R *5 . D+ 0 
THT0=R*1. 30+1 
X (N2, 1 ) = 0. D + 0 

X (N2 , 2 )=NING*X (Ml, 1 ) +NTN0*X ( N 1 , 2 ) -F I VC*X ( N1 , 3 ) +R* X ( N 1 H 

N M 3 = N - 3 

DO 1 K = 1 , NM 3 

KP1=K+1 

KP2=K+2 

KP 3 = K + 3 

ADD=THTO* ( X(N1,KP1)+X(N1,KP2))-R*(X(N1,K)+X(N1,KP3J) 

1 X (N2, KP2)=X (N2 ,KP1 ) +A00 

X (N2,N )=X(N2,N-1)+NIN0*X( N 1 , N ) + NTNO *X < N1 , N- 1 ) -FIVC*X(N1 
1+R*X (Nl,N-3) 

RETURN 

ENC 
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Appendix C 



Examples 

The listing on the next page is of a MAIN program which 
calls RODISK twice to solve two different problems. On an IBM 
360/67 compile time (for the MAIN, RODISK, and the ancillary 
subroutines) was 12 s, link time was 2 s, and execution time 
for both broblems was 1.5 s. In both problems we used N = 4l 
and iterated 10 times. 

The first problem is of type 2 (solid disk) with b = 10, 
t = l/(1.6+.008r 2 ) , v = 0.3, yea 2 = 120, a r (b) = 14000, and 
EotT = 25105 + 1300 log g (t) - r 2 (233+l6t) Units are inches and 
pounds. The problem was made up from the exact solution 
a r = 9000 + 50r 2 ; a 0 = a r + r 2 (120+l6t) 

The results, shown on page 18, show evaluations for the 
stress components which are correct within 0.03 psi even though 
convergence was complete to only about 4 digits as indicated by 
the sequence of values above the final tabulation; these are 

values of a (0 )/o (b) and of B. 

r r 

The second problem is of type 1 with a = 2.57, b = 5.15, 
a (a) = 18205, a (b) = 22000, t = 0.1493 r"^ 2 , T = 60 - 1.6r 2 , 
v = 0.3, and Ea = 194.3. The units are inches and pounds. 

Since the thickness variation is a power law, the theoretical 
solution given in the body hereof may be used to obtain the 
theoretically exact solution 

a = -113-95r 2 + 15832^ - 11708r q 
r 

a n = 122.80r 2 + 13801^ + l4310r q 
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C J . E . BROCK 1 MAY 1976 

C THIS IS A MAIN PROGRAM TO TEST MY SUBROUTINE RODISK 

IMPLICIT REAL*8(A-H,0-Z) 

REAL*S X(2J , 131 ) 

INTEGER K P ( 4 ) 

COMMON X,N 

COMMON /ONE /ARAO, BRAD ,SRA , SRB , T EEBE E , POI S , KP 
C THE FIRST PROBLEM IS FOR A SOLID DISK 

KP ( 1) =2 
KP (2 ) =4 
KP (3) =10 
K P ( 4 ) =1 
BRAD= 1 . D+l 
SR 8 = 1.40+ 4 
N = 41 
E NM=N- 1 
DO 20 1=1 ,N 
c I M = I - 1 

R =EI M*8RAD/ENM 

TH IC K = 1 . D+0 / ( 1 . 6D+0+8.D-3*R**2) 

X (2) , I )=THICK 

X ( 2,1 ) =12 0. D+0+BRAD**2/SR9 

W = 2. 510 50 + 4+1. 30+ 3"DL0G(THICK)-(R**2) *(2 . 33 0 + 2+ 1 . 60 + 1*TH I CK 
X (3, I ) =W/SRB 
20 CCNTINUE 
S = X ( 2 C , N ) 

TEEBEE=S 

CALL CIVS(20,1,S) 

CALL ROOISK 

C THE SECOND PROBLEM IS FOR AN ANNULAR DISK 

TcEBEE=l. 4930-1*5. 150+0** (-4.20-1) 

ENM=N -1 
KP (1 ) =1 
ARA0 = 2. 570+0 
3R A0= 5 . 15 D + 0 
SR A= l . 82050+4 
SRB = 2 . 20 + 4 
BETA=5. 32473471C-1 
00 30 I =1 ,N 
2 1 M= I - 1 
Y=EI M/ENM 
X ( 5 , I )=Y 
X (6, I )=Y 

R =ARAD+( BRAD- A RAD ) *Y 
X ( 4, I ) =R/BRAD 

THICK=1.4930-1*R**(-4.2D-1) 

X( 1,1 ) =7H I CK/TEEBEE 
X ( 2 , I ) = BETA 

X(3,I)=(6. 0+1-1 .6C+0*R**2 )*1 .9430+2 /SRB 
30 CONTINUE 
CALL RODISK 
S T CP 
END 
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where p = .29171 and q = -1.87171. The computer results agree 
with five digit accuracy even though the sequence of iterant s 
(A,B) shows only four digit convergence. 

The same two problems were also worked with N = 101 and 14 
iterations. The execution time went from 1.5 s to 3*0 s. The 
maximum change in any stress value was 0.3 psi. From these and 
other problems it may be concluded that there are no difficulties 
of accuracy, computer storage, or execution time. 
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